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ABSTRACT 


Euler codes have been developed to study axisymmetric and 
two-dimensional internal flows . A time-dependent approach has been 
used to predict steady state flowfields. The steady state solution is 
achieved asymptotically. The method is based on backward finite 
differencing in time. A least-squares finite element scheme for first 
order systems of partial differential equations in space is applied to the 
Euler equations. The scheme minimizes the L 2 -norm of the residual within 
each time step. The method naturally generates numerical dissipation 
proportional to the time-step. An implicit method employing linear 
elements has been implemented and has proved robust. Three separate 
codes have been developed and validated for 2-D case in non- 
conservative variables, 2-D conservative variables and axisymmetric 
case. The various cases that have been tried include shock reflection from 
a wall, supersonic flow over a ramp, two-dimensional nozzle with 
supersonic inlet flow, 2-D mixed flow nozzle with and without shocks, 
axisymmetric no zzl es and flow through a turbine cascade. Flowfields in 
some of the cases involved both the subsonic and supersonic regimes. The 
shock capturing capability of the codes has also been demonstrated. The 
results have been compared wherever possible with experimental or 



computational results from various sources and good correspondence has 
been obtained. The method is computationally expensive especially when 
shocks are involved or when the grid is relatively fine. Matrix solvers like 
Gauss elimination, Cholesky’s LDL T factorization method as well as 
frontal solver were used. Frontal solver reduced the RAM requirements 
but increased disk space requirement in comparison to the Cholesky’s 
method but there was no appreciable time saving observed between the 
two for smaller systems. 



CHAPTER 1 


INTRODUCTION 

Propulsion problems encountered in the real life are complex 
due to interplay of heat and mass transfer, fluid mechanics, chemical 
kinetics, thermodynamics and turbulence. Such problems are not tractable 
analytically. It is in this context that numerical methods have assumed 
eminence. Potential of numerical methods could be exploited when 
algorithms are developed and chosen judiciously as well as adequately and 
fast computing system is available. 

The internal flows encountered in different components of a 
propulsion system such as intakes, nozzles and turbomachinery cascade 
flows are frequently characterized by mixed subsonic and supersonic 
regimes and at times presence of shocks. Such flow fields cannot be 
studied by either method of characteristics or finite difference space 
marching technique due to mixed nature of the flow. Time marching 
technique has been used for quite some time now to simulate such 
flowfields because it renders the equations hyperbolic with respect to time 
and so for the entire domain of interest the nature of equations remains the 
same. The pseudo-transient is used to obtain the steady state results. 

The Navier-Stokes equations govern the flows commonly 
encountered in both internal and external applications. Computing a 



solution of the Navier-Stokes equations is frequently impossible or atleast 
impractical and in many of these applications, unnecessary. For sufficiently 
large Reynolds numbers the important viscous effects are confined to a 
thin boundary layer near the surface of a solid boundary. As a 
consequence , the inviscid (non-viscous, non-conducting) portion of the 
flowfield can be solved independently of the boundary layer. Ofcourse this 
is only true if the boundary layer is very thin compared to the characteristic 
length of the flowfield so that the interaction between the boundary layer 
and the inviscid portion of the flowfield is negligible. The reduced set of 
equations valid only in the inviscid portion of the flowfield, obtained by 
dropping both the viscous terms and the heat transfer tenns from the 
complete Navier-Stokes equations can be numerically solved using much 
less computer time than is required for the complete Navier-Stokes 
equations. These equations are referred to as the Euler equations. In the 
present investigation a few typical propulsion problems involving internal 
flows have been studied using time dependent finite element method. Euler 
equations have been solved numerically based on the least-squares 


formulation. 



1.1> LITERATURE SURVEY 


A> FINITE DIFFERENCE AND FINITE VOLUME TIME 
MARCHING TECHNIQUES 

In a mixed flow flow field the nature of equations changes from 
elliptic to hyperbolic or vice-versa and this creates difficulties in the 
numerical solutions. Hence time dependent method is frequently used to 
analyze such flow fields. The time dependent method uses a pseudo- 
transient approach whereby the steady state is obtained asymptotically. 
Using this approach the equations remain hyperbolic with respect to time 
for both subsonic as well as supersonic regimes and so the difficulties 
faced previously are avoided. Very efficient shock capturing schemes have 
been developed for these equations for internal and external flows. The 
discontinuities are obtained as part of the solution. The advantages of time 
marching are: 

1> The same code works for the all flow regimes whether subsonic or 
hypersonic, viscous or inviscid; 

2> The technique can be easily implemented and modified for use on 
parallel computers; 

3> The same computer program can be used to analyze internal as well as 
external flows by using the appropriate boundary conditions; and 



4 > For complex flow geometries and flow fields the technique is robust. 

The time marching schemes can be either explicit and implicit. The 
Lax -Wendroff scheme and its variants were the first explicit schemes to 
be developed. After this came the predictor-corrector scheme of 
MacCormack [1] . Jameson developed and implemented the Runge-Kutta 
multistage methods for Euler equations in 1981 [1]. This method is 
especially attractive when implemented with implicit residual averaging, 
local time-stepping and multi-grid finite volume formulation. In case of 
implicit schemes the time-step is not restricted by the CFL limit. A 
comparison between the explicit and implicit schemes for turbomachinery 
flows by Hall [1] showed that the Runge-Kutta method gives slightly 
better performance as the artificial viscosity is introduced and controlled 
by the user. 

B.> FINITE ELEMENT METHODS 

One of the major problems faced in finite difference and finite 
volume methods in CFD is the generation of grids for three dimensional 
geometries realistic geometries of aerospace vehicles. The generation of a 
single grid that discretizes the entire flow region for a complex region is a 
very complicated task. Even for axi-symmetric geometries structured grid 
generation for multiple connected domains such as that of an airbreathing 



rocket, often leads to problems such as improper resolution in some 
regions and excessive grids in some other regions. These difficulties have 
led to recent research into alternate approaches for handling complex 
geometries such as zonal approach, unstructured finite volume and 
unstructured finite element approaches. 

Finite element methods are rather new to CFD and hence immature 
for flowfield analysis in comparison to finite difference method. The finite 
element method was originally developed by aircraft structural engineers 
in 1950's to analyze large systems of structural elements in the aircrafts. 
Turner, Clough, Martin and Topp presented the first paper on the 
subject, followed by Clough and Argyris , among others [2], Applications 
of finite element method to non-structural problems such as fluid flows 
was initiated by Zienkiewicz [2]. 

Today various theories of fluid behavior are available which 
encompass virtually any type of phenomena of much immediate practical 
interest. However, there remains a surprising number of unsolved 
important problems in fluid dynamics due to difficulties encountered in 
most of the conventional analytical and numerical methods. In fluid 
dynamics, a choice of Eulerian coordinates renders the resulting 
governing equations non-linear, in general, and thus analytically difficult 
to solve. The most widely used numerical method of overcoming these 



difficulties has been the method of finite differences. Among other 
popular methods are the variational methods and methods of 
weighted residuals. Unfortunately, variational principles often cannot be 
found in some engineering problems, particularly when the differential 
equations are not self-adjoint. Weighted residuals are applied in the 
methods of Galerkin, Least-squares, and Collocation. In case of least- 
squares method, the physical behavior of the system is generally 
described by linear or lower order functions. 

The success of upwind finite differencing and related artificial 
dissipation methods motivated studies of analogous upwind finite 
element methods; similarly, the idea of the Lax-Wendroff scheme in 
finite differencing has motivated studies of Taylor-Galerkin finite element 
algorithms. For example, the streamline upwind Petrov-Galerkin 
method [3], the Taylor-Galerkin method [4-6], the Taylor-Galerkin 
method with flux-corrected transport (FCT) [7-9], block relaxation via 
Godunov's method [10] and the characteristic Galerkin method [11] have 
been developed and applied with considerable success to these problems. 
Recently attempts are being made on the parallelisation of Finite Element 
programs to increase the overall efficiency of the computation [12]. 

In the least-squares method, generally linear order elements 
are used and they give reasonable results. Moreover, the stiffness matrix 



produced for least-squares formulation is symmetric and positive definite 
so matrix solvers like the Cholesky’s method can be used resulting in the 
program being more efficient. Jiang and Carey [13] have used the time- 
dependent version of this method on Euler equations and demonstrated its 
robustness in shock capturing. 

1.2> PRESENT INVESTIGATION 

In the present investigation, three separate codes have been 
developed based on Ref. [13], These codes are modular in nature and 
hence easily alterable as and when necessary. They are capable of 
simulating both external as well as internal flows, depending on the 
boundary conditions. The codes have been developed for: 

a. ) 2D case in non-conservative variables; 

b. ) 2D case in conservative variables; and 

c. ) Axisymmetric case in conservative variables. 

The code can only simulate inviscid flow fields, being based on 
Euler equations but can be easily extended to study viscous flow regimes 
as well. It uses the time-dependent least-squares finite element 


formulation. 



In solving the matrices produced by the formulation Gauss 
Elimination, Cholesky’s LDL T factorization as well as Frontal solver were 
used but generally Cholesky’s method was used. 

A number of steady state cases have been studied. The different 
problems studied are: 

1. ) shock reflection problem based on Ref. [13] using the 2D 

nonconservative variables code; 

2. ) Supersonic flow over a 20 degree ramp based on Ref.[13] 

using 2D non-conservative variables code; 

3. ) B3 nozzle for pressure ratios of 9.0 & 1 .98 [1] using the 2D 

conservative variables code ; 

4. ) Axisymmetric nozzles [Ref. 14,15] without shock using the 

axisyrnrnetric code; 

5. ) 2D nozzle with a half angle of 15 degree using the 2D 

conservative variables code; and 

6. ) Turbomachineiy cascade flow [Ref. 16] using the 2D 


nonconservative variables code . 



CHAPTER 2 


METHODOLOGY 


1 . Governing Equations 

The compressible Euler equations, which form the basis of this 
thesis, can be written in conservative form as follows: 

™- + — -I G " L ) j0 JL_ (2.1) 

dt dx y 10 dy y >0 


where j 0 = 0 or 1 for either two-dimensional or axisymmetric flow. 
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The conservation form ( or conservative variables ) implies that the partial 


differential equations (PDE) having this form have the property that the 
coefficients of the derivative terms are either constant or, if variable, their 


derivatives appear nowhere in the equation. For a PDE which represents a 
physical conservation statement, this means that the divergence of a 
physical quantity can be identified in the equation. In Navier- Stokes 



equations when the dependent variables are in the form: 
p , pu , pv , pe 

then the equation is said to be in conservative form. On the other hand, if 
the dependent variables are: p , u , v , e 

then the equation is said to be in non-conservative fonn. The finite 
representations of the PDE should be such that a good approximation of 
the conservation laws implied by the PDE is maintained. Almost all the 
finite element analogs can be generally assured of the conservative 
property because of the veiy nature of the control volume built into the 
integral formulation in FEM. 

2> Technique 

The steady supersonic flowfields are governed by 
hyperbolic differential equations, whereas steady subsonic flowfields 
are described by elliptic differential equations. 

There are various cases where the internal and/or external flow 
involves both the subsonic and supersonic regimes. Flows through 
converging-diverging nozzle and over a blunt body may involve both the 
regimes. The fact that the nature of the governing equations changes 
from elliptic to hyperbolic across the sonic line causes severe 
mathematical and numerical difficulties at such a level that steady flow 


in 



solution of the subsonic and supersonic regimes are usually treated 
separately and differently. The time marching technique for the solution of 
steady flows was developed in 1966 by Moretti and Abbett [20], which 
was published as the first practical solution for the supersonic blunt 
body problem. Since then time-marching also called as time dependent 
solutions have become an important segment of computational fluid 
dynamics. 

The unsteady equations of motion are Hyperbolic with respect to 
time, regardless of the fact whether the flow is locally subsonic or 
supersonic. Hence the flowfields lend themselves to a well posed 
initial value problem with respect to time. Therefore, the time marching 
technique becomes very powerful tool for the solutions of such mixed 
flows, being uniformly valid throughout the flowfield. 

In flowfields involving shock waves, there are sharp, 
discontinuous changes in the primitive flowfield variables such as 
pressure, density, velocity, temperature etc. across the shocks. Many 
computations of the flows with shocks are designed to have the shock 
waves appear naturally within the computational space as a direct result 
of the overall flowfield solution, i.e. as a direct result of the general 
algorithm, without any special treatment to take care of shocks 
themselves. Such approaches are called shock capturing methods. This 



is in contrast to the alternate approach, where shock waves are 
explicitly introduced into the flowfield solution, the exact Rankine- 
Hugonoit relations for changes across shock are used to relate the flow 
immediately ahead of and behind the shock, and the governing flow 
equations are used to calculate the remainder of the flowfields. This 
approach is called the shock fitting method. 

The shock fitting teclmique often breaks down when shocks 
develop spontaneously within the fluid. 

Here a time dependent finite element method with shock capturing 
has been used. Finite element methods are well suited to practical 
problems posed on the complicated domains and are now extensively 
applied in solid mechanics, heat transfer, fluid flow etc. The earlier 
difficulties encountered in extending the methodology beyond elliptic 
boundary value problems and to convection dominated flows have 
been substantially resolved. Here we begin by considering the backward 
implicit time differenced formulation. 



3. 


Mathematical formulation 


The three different cases of two-dimensional non-conservative 
variables, two-dimensional conservative variables and axisymmetric 
conservative variables are treated separately as follows: 


a> The two-dimensional non-conservative casef 131: The two- 

dimensional unsteady compressible Euler equations in non-conservative 
form can be represented as a first-order system as: 


dU 

dt 


+ A, 
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+ A , = 0 

2 dy 


in Q x (0,T), 


MU = g in T g x (0,T), (2.3) 

U = U 0 in Q for t = 0 

where U T = (p u v p), M is a boundary operator, g is a given vector- 
valued function, T g is that part of the boundary T where essential 
boundary conditions are applied, Q is the spatial domain of the problem, t 
is time and 
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Here p is the density, (u,v) are the velocity components, p is 
the pressure and y is the specific heat ratio. 

For a given time step At = t I1+1 - f we linearize the problem 
by setting Ai=A](U 11 ) , A 2 = A 2 (U n ). Backward differencing leads to the 
implicit time-differenced problem 


B+1 TT " ' A * Am dU " 1 ■+ AtA” — — = 0 
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dx 
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(2.5) 


Introducing AU = U n+1 - U 11 we rewrite the above equation as 
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= 0 ( 2 . 6 ) 


The basic least-squares method for the system (2.6) amounts to 
minimizing the L 2 - norm of the residual R for admissible AU in (2.6); i.e. 
minimizing the objective functional 


O 0 = J R T Rdxdy 


(2.7) 


with 
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( 2 . 8 ) 



Taking variations with respect to AU and setting the test 
function V = 8AU, § <t> 0 = 0 leads to the least-squares weak statement: find 
AU e S = = 0 on where H^D) denotes the usual 

Hilbert space, such that 
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Next we introduce a finite element discretization. Let N e be 
the number of nodes for an element and ¥ j denote the element basis 
functions. The approximation and the test function are 



( 2 . 10 ) 


V h (x,y) = v F l (x,y)E (2.11) 

where A(p,u,v,p)/ are the nodal values at the j* node and E 
is the 4x4 identity matrix. Substituting AUh, Vh for AU, V in (2.9) and 
evaluating the integrals, we have 

K(AU) = F (2.12) 

where AU is the global nodal vector. The global matrix K and the global 
vector F are assembled from the following submatrices and subvectors 
respectively: 



( 2 . 13 ) 


K l = \(Lq>,) T lAPjdxdy 
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where L<p i = <p t E + A/<z o ix A" + At<p l v A% 


It can be noted that the stiffness matrix K is symmetric, 
positive and definite. This is an important feature of the least-squares 
method. Only the lower half matrix need be generated and stored if 
Cholesky LDL T factorization is used. 


b.> Two-Dimensional conservative formfl31 : The Euler 

equations governing two-dimensional compressible inviscid flows can be 
written in conservative form as: 
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Q - Q 0 in Q for t = 0 
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in which e is the total energy and for the case of a perfect gas the equation 
of state is 
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Now we convert this to the following modified form: 
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(2.18) 


It has to be noted that the A-matrices (both 1 & 2) in this case 
are different from those of the non-conservative case. The new A-matrices 


are given below: 
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Using the same procedure as before, we construct a 
corresponding least-squares weak statement similar to (27) where the 
matrix operators in (23), (24) are again evaluated at the previous time level. 
This implies that an approximate conservative fonn can be reconstructed. 
We have tire a conservative least-squares weak statement: 

find A QeS= {(//’ (n)) 4 ; M(AQ) < 0 on T g J such that 
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here 


dF 


dx 


A „ dQ” dG ■ „ dQ 

1 dx ’ dy ^ dy 


The unknown here is AQ, which is the increment in 
conservative variables for time step At. We may use the conservative 
variables at the previous time step to calculate the nodal values of 
components of flux F n and G n . 


c.> The Conservative Axisymmetric case: For this case the 


equations are as follows: 
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Here again almost the same analysis as that of the case (b) is 
repeated. Pressure P is represented by: 


P = (r-^\p pu p\> 





(2.23) 
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The above equation can be reduced to the following form: 


When the analysis is earned on the above equation the first three terms are 
treated as in the 2-D conservative variable case. Only the third new term 
has also to be taken into account. The A\ and A 2 matrices are the same as 
in the case (b). Here the residual R is 
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The M and S matrices are as follows: 
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Now again minimizing the L 2 -norm we get the corresponding weak form 
and thereby the K and F matrices. 


CHAPTER 3 


A NOTE ON THE COMPUTER PROGRAM 

The mathematical formulation has already been discussed in 
the last chapter. Utilizing that formulation three separate codes have been 
developed. These codes are modular in nature and as such with the 
exception of a few subroutines the codes are same for all three cases. 
Given below is a listing of the subroutines along with the function of each 
for the convenience of any prospective user. In all there are 16 
subroutines: 

1> NODEDT : In this subroutine the geometry of the domain is 
specified and it generates the nodes and the connectivity matrix. Here the 
number of nodes and the number of elements is supplied. 

2 > INITDT : In this subroutine the initial values of the dependent 
variables are specified and the number and location of the gauss points has 
to be given for numerical integration. The weights for different orders of 
integration at different gauss points are also given. 

3 > WAJAZ : This subroutine is called from the subroutine 
ELEMENT. It supplies a component of the elemental stiffness matrix at a 
particular gauss point. It has to be supplied the number of element being 
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computed and the gauss point for which the value is being presently 
calculated. 

4 > SABAZ : This subroutine is also called from the ELEMENT. 
This requires the element number and the gauss point and gives the value 
of the elemental force vector at this location. 

5> MULT : This subroutine is simply for multiplication of different 

matrices required in the computation of elemental force and stiffness 
vectors. 

6 > TRANS : This subroutine simply gives the transpose of the 
matrix supplied to it. 

7> ERROR : This subroutine is called from different subroutines if 
there is an error. It prints an error message. For example, when the 
jacobian becomes equal to zero (which it should not be) or when the 
global stiffness matrix which should be symmetric but is not due to some 
error. 

8> GKSS : This subroutine is one of the most important one’s. It 

generates the elemental force and stiffness matrices by calling ELEMENT 
and then assembles them to yield the global stiffness and force matrices. 

9> SHAPE: This subroutine supplies the shape functions for the 
four noded quadrilaterals that have been used in this analysis. 



]0.> JACOBI : This supplies the Jacobian for the isoparametric 
transformation of the elements and also checks that the jacobian is not zero 
otherwise it would imply that there has been faulty grid generation. 

1 1 .> ELEMENT : This subroutine has already been mentioned above and 
it is called from GKSS. It requires the element number. It does the 
numerical integration and supplies the elemental stiffness and force 
matrices. 

12> APLYBC : As is obvious from its name, this subroutine applies the 
boundary conditions. 

13.> CALFUN : This subroutine calls the matrix solver like the gauss 
elimination or the Cholesky’s factorization, as the case maybe. 

14> PRINTS : This one prints the results. All the computed variables 
are supplied to this subroutine through the common block and whichever is 
desired can be printed. 

15> CHOLDC : This is a part of the two subroutines that do the 
Cholesky’s factorization and then solve the matrix. This does the 
decomposition part. 

16> CHOLSH: This solves the matrix to yield the results of the 
algebraic equations that result from any numerical method. 

The time step, maximum number of iterations as well as the 
tolerance factor are specified in the main program. For different cases (i.e. 



whether 2D or axisymmetric, conservative or non-conservative) the 
subroutines SABA2 and WAJAZ have to be replaced with the appropriate 
ones. For different geometries the NODEDT will have to be modified. If 
frontal solver is used then GKSS has to be deleted because frontal solver 


does the assembly itself and in a different manner. 



CHAPTER 4 


RESULTS & DISCUSSION 

Seven cases have been run using the developed codes for 
different boundary conditions. They are discussed below in serial order: 
a.> Shock Reflection Problem : The isobars for the standard test 

problem corresponding to the reflection of a shock from a wall are 
depicted in Fig. 1. On the upper boundary of the flow domain p=1.7, 
u=2.6185, v=-0.5082, p=1.5282 and on the upstream boundary p=1.0, 
u=2.9, v=0.0, p=0.7143, so a shock emanates from the upper left comer. 
This shock is reflected at the lower wall where v=0 (symmetry condition is 
applied here) and the downstream boundary conditions remain free for 
outflow. The initial data were prescribed as constant at values given on the 
upper boundary and the specific heat ratio is y=T .4. In the calculation a 
uniform 60x20 mesh of bilinear elements was used. The solution was 
integrated with a time step At=0. 33333 until an essentially steady state 
was obtained in 26 timesteps. The results were obtained using the 2-D 
code with non-conservative variables. The isobars indicate that the shock 
captured is somewhat smeared but the shock angles and the pressure 
values at the exit (i.e. after the shock reflection has already taken place) 



show excellent correspondence with the theoretical results. It can be 
inferred that the flow modelling has been correctly done. This problem 
was taken from Ref. 13. 

b.> Supersonic Flow over a Ramp : The second case studied has 

again been taken from Ref. [13]. Here a supersonic flow with a Mach 
Number of 3, over a 20° ramp has been studied using the 2-D non- 
conservative variables code. The gas enters with uniform flow conditions 
through the left boundary of the domain and an oblique shock develops at 
the root of the ramp. The inlet parameters are u=3.0,v=0.0,p=1.0 and 
p=0.7143 and the same boundary condition was applied on the upper 
boundary. The mesh contains 824 bilinear elements and the computed 
pressure contours are illustrated in Fig. 2. In the calculation the initial data 
were prescribed as constant at the value given on the left boundary and the 
time step was At=0. 33333. The steady state was obtained in 30 time steps. 
The result obtained here (for eg. the shock angle and the pressure at the 
exit ) are again in close agreement with the theoretical results, although the 

shock, as in the previous case, is somewhat smeared. 
c > m Nr> 7 ? 1 fv Fig.3 shows the dimensions of B3 nozzle. This 

geometry, input data and the experimental data (for comparison) were 
taken from Ref. [ 1 ] . This nozzle was studied for two pressure ratios of 9.0 



& 1 .98. The velocity components(u&v) and the density are specified at the 
inlet(p=1.0,u=0.21,v=0.0) and the pressure ratio at the outlet(Po/p=9.0). 
For the first pressure ratio of 9.0 the flow field has no shock. The analysis 
was done using 30x30 bilinear elements. The results were obtained using 
the 2-D conservative variables code in 1 5 time steps and the time step was 
the same as in previous case. Fig. 4 shows the Iso-Mach contours and Fig. 
5 shows the iso-bars for this case. Figures 6 & 7 show' the centre-line and 
wall pressure distribution respectively. They are in good agreement with 
the experimental data marked on the figures. Figures 8 & 9 show the 
centre-line & wall Mach number distributions respectively. The Mach 
number in actual case must be zero at the wall but since an inviscid code 
has been used so it is not. 

Figures 10 & 11 show the Iso-Mach & iso-bars for the same 
nozzle but for the pressure ratio of 1.98. At the inlet p=1.0,u=0.21,v=0.0 
and at the outlet Po/p=1.98. The time step is again At=0. 33333 and the 
results were obtained using the 2 -D code in conservative variables. The 
shock that exists for this pressure ratio at about x=1.2 may not be clear 
from the above figures but can be clearly observed from the centre-line & 
wall distributions of Mach number and pressure. For this case figures 12 & 
13 show the centreline & wall pressure distributions respectively. In 



these figures (12 & 13) the experimental data was available so has been 
shown on the plots and the two results (computationally obtained and the 
experimentally obtained ones) indicate a good matching. No results were 
available for the mach number distributions shown in figures 14 & 15 
which illustrate centre-line & wall Mach numbers respectively. 
d> Axisymmetric Nozzle : Fig. 16 shows the nozzle geometry taken 
from Ref. [15], 0i=29° , Y B = 1.2, e B =30°,X B =3.6, Y E =3.7 and X E =9.4. 
This case was run with the pressure ratio (p/P 0 )of 0.01 and the input 
uniform flow of Mach number of 0.1. On the wall the flow tangency 
condition was applied and on the lower boundary the symmetiy condition 
was used. The time step was the same as in all previous cases and 40x10 
mesh was used. The results were obtained in 24 time steps. Fig. 17 & 18 
show the iso-machs and iso-bars respectively. Wall pressure distribution is 
shown in Fig. 19. Fig. 20 shows the centre-line pressure distribution. The 
results from Ref. [15] obtained using a finite volume discretization are also 
shown alongside and the results are once again in good agreement. In 
figures 21 & 22 the wall & centre-line Mach. no. distributions are given 
and they also match properly with the results of Ref. [15]. 
e> 2-D Nozzle : Fig. 23 shows the mesh and the geometry used 

for this nozzle. This was a purely supersonic case with no transition from 
subsonic to supersonic domain. The mesh was very coarse with only 70 



elements which were again bilinear. The half angle was 1 5°. and the case 
was run for the inlet pressure (non-dimensional) of 0.71 and the 
corresponding Mach number with a uniform inlet flow. The results match 
with one-dimensional results. Fig. 24 shows the wall pressure distribution 
for this case. Fig. 25 shows the Mach number distributions for different 
axial stations. The time step used was the same as in previous cases and 
the solution was obtained in 8 time steps. 

£> TND-2579 Nozzle :The data for this was taken from Ref. [14]. The 

geometry and the mesh used for this nozzle is shown in Fig.26. This again 
is a purely supersonic case and the axisymmetric code was used to obtain 
the results. The specific heat ratio for this case was y=1.2. The inlet 
uniform flow had a Mach. no. of 1.05(p=1.0,u=1.05,v=0.0). Here a mesh 
of 35x8 was used. The time step was At=0.5 and the solution was obtained 
in 20 time steps. Fig. 27 & 28 show the centre-line & wall pressure 
distributions for this case. Fig. 29 shows the centre-line Mach. no. 
distribution for this nozzle. Fig. 30 shows the wall Mach. no. distribution 
and for this figure we have method of characteristic results [14] to 
compare. Once again the comparison is excellent. 

g> Cascade Flow : For this case a mesh of 40x8 was used. The 

airfoil used has been shown in Fig.31. A typical mesh used in the analysis 



of gas turbine cascades has been shown in Fig. 32 (this is not the actual 
mesh). The time step was At— 0.33333 and the steady state was reached in 
16 steps. The inlet angle was 44.5° and the outlet flow angle was 21.8°. 
The cascade pitch was 0.825. The inlet flow Mach, number was 0.32 and 
the pressure ratio imposed at the exit was the one corresponding to the 
Mach number of 0.58. A periodic boundary condition was imposed on the 
portions of the computational domain such as: AB & IH and DE & GF 
(this implies that the values of the computed variables at the lower 
boundaries IH & GF are updated by the values at the upper boundaries AB 
& DE respectively after each iteration). Fig.33 shows the centre-line 
Mach no. distribution and the results show a good matching with those 
obtained by the streamline curvature analysis. Fig. 34 shows the pressure 
coefficient for this case and they appear to be slightly off compared to the 
results of Ref. [16]. The results from the present analysis depart from those 
of the Ref. [16] especially near the leading and trailing edges. This may be 
due to the fact that in the reference quoted above the results have been 
arrived at by repeating the computation around the edges after an initial 
inviscid analysis followed by a transitional profile boundary layer and 
wake analysis. 

* 

The results have demonstrated the robustness of the codes. The 
advantage of the least-squares formulation is that reasonable results can be 



obtained by lower order elements. In this formulation artificial viscosity 
appears naturally and the matrices are symmetric (thereby improving the 
efficiency of the computation). But alongside the shortcoming in this 
formulation (L 2 ) is that it is only first order accurate in time and if we try to 
improve the accuracy to second order then this may lead to oscillations 
(since we are investigating steady state cases so time accuracy is not 
relevant). For higher order elements this formulation may result in the 
appearance of non-physical oscillations in the solution. 



CHAPTER 5 


CONCLUSIONS AND SCOPE OF FURTHER WORK 

5.1> CONCLUSIONS: In order to predict the internal flowfields 

Euler codes have been developed based on least-squares finite element 
method. A time-dependent approach has been used to overcome the 
problems encountered in studying, computationally, flows involving both 
subsonic and supersonic regimes. Iso-parametric (bilinear quadrilaterals) 
elements have been used. The codes were developed for two-dimensional 
as well as axisymmetric cases. For two-dimensional case , two separate 
codes were developed, one each for non-conservative variables and 
conservative variables. The axisymmetric code uses conservative 
variables. 

The 2-D non-conservative variables code was first used on 
the shock reflection problem. Proper shack capturing took place and the 
computed results matched well with the theoretical shock angles and the 
after shock pressure. 

The non-conservative 2-D code was also used to study 
supersonic flow over a ramp. Once again the results showed good 
correspondence with the theoretical results. 



The 2-D conservative variables code was used on the 
B3 nozzle for two different conditions, one invoking shock and the other 
without it. The results were slightly off but trends were the same in both 
the cases. 

The results for the axisymmetric nozzle, without shock were also in 
agreement with those of a finite volume analysis [15]. This case also 
involved both subsonic and supersonic flow speeds. 

The axisymmetric code was also run for NASA TND 2579 
nozzle. This was a purely supersonic case and the results were matched 
with method of characteristics results and once again there is veiy good 
correspondence. 

The 2-D conservative code was used on a two-dimensional 
nozzle with supersonic inlet conditions. The results were compared with 
one dimensional analysis and the results were almost the same as the 
theoretical results. 

A subsonic flow through a gas turbine cascade was analyzed 
using the non-conservative 2-D code and the results were slightly off 
compared to those of streamline curvature analysis. 

The results proved the capability of the codes to predict the 


flowfields for internal flows, with and without shocks. 



5.2.> SCOPE OF FUTURE WORK: The entire analysis so 

far has been done with bilinear quadrilaterals. The same analysis can be 
repeated with higher order elements. Same 2-D case can be run for both 
conservative and non-conservative variables mid the results compared. 
Because of the paucity of time the number of cases tried were limited. 
More complex cases (for e.g., mixed flow through the turbine cascades) 
can be tried. Different new geometries can be tried, especially those for 
which the grid generation is more complicated for e.g. internal 
compression intakes and hypersonic nozzles with flaps. The use of 
adaptive grid generation can be made to observe how the quality of results 
alters. New subroutines can be added to include the viscous effects in the 


program. 
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Fig, 1 . Iso-bars for Shock Reflection problem 
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Fig. 3. 2D convergent-divergent (B3) nozzle geometry 


24.6 








X/Le — -> 

Fig. 6. Centre-line Pressure distribution for B3 nozzle (Po/p=9.0) 
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Fig. 7, Wall Pressure distribution for B3 nozzle (Po/p=9.0) 
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Fig. 8. Centre-line Mach no. distribution for B3 nozzle (iyp=9.0) 
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Fig. 9. Wall Mach No. distribution for B3 nozzle (iyp=9.0) 




Fig. 10. Iso-Machs for B3 nozzle (iyp=1.98) 


49 




Fig. II. Iso-bar? for B3 nozzle (IVp=1.98) 



50 



CENTRE-LINE PRESSURE 
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Fig. 13. Wall Pressure distribution for B3 nozzle (Po/p=1.98) 
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CENTRE-LINE MACH No. 
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Fig. 16 . Axisymmetric nozzle geometry 
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Fig. 17. Iso-machs for axisymmetric nozzle 
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Fig. 18 . Iso-bars for axisymmetric nozzle 
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Fig. 21 . Wall Mach No. distribution for Axisymmetric nozzle 
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Fie. 25. Mach No. distribution at various axial stations for 2D nozzle 
(half angle 15°) 
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Fig, 26, TND-2579 nozzle geometry (Axisymmetric) 
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Fig. 29. Centre-line Mach No. distribution for TND-2579 nozzle 
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Fig. 30. Wall Mach No. distribution for TND-2579 nozzle 
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Fig. 31. Cascade flow Airfoil geometry 
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